# setup
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
np.random.seed(23)

# setup: has to be in a separate cell for some reason
plt.rcParams['figure.figsize'] = [10, 5]

Individual-based models

Let's start with a simple model of the size of population of organisms. We'll ignore sex (so, suppose these are hermaphroditic).

Population growth

Each year, each individual will

1. maybe reproduce, and then
2. maybe die.

To make a model we need to decide exactly how those steps happen. How about:

1. Each individual produces a random, Poisson number of offspring, with mean $\lambda$, and then
2. survives with probability $p$ (or else, dies).
3. Every individual does this independently of everyone else.

def exp_gen(N, lam, p):
    """
    N: current number of individuals
    lam: mean number of offspring per capita
    p: probability of survival per capita
    Returns: the number of individuals in the next time step
    """
    births = np.random.poisson(lam, N)
    survivors = (np.random.uniform(size=N) < p)
    return sum(births) + sum(survivors)

# run five simulations
N = np.full((30,5), -1, dtype='int')
N[0, :] = 1
for t in range(1, len(N)):
    for k in range(N.shape[1]):
        N[t, k] = exp_gen(N[t-1, k], 
                          lam=0.5, # mean offspring / year
                          p=0.75) # mean lifetime is 1/(1-p)

assert(np.min(N) >= 0)

fig, ax = plt.subplots()
ax.plot(N)
ax.set_xlabel("time")
ax.set_ylabel("pop size")
ax.set_title("Exponential growth")
# ax.set_yscale("log")

Text(0.5, 1.0, 'Exponential growth')

N

array([[   1,    1,    1,    1,    1],
       [   2,    1,    1,    1,    1],
       [   2,    2,    0,    1,    1],
       [   1,    3,    0,    2,    1],
       [   2,    3,    0,    4,    3],
       [   1,    4,    0,    7,    4],
       [   0,    6,    0,   12,    6],
       [   0,    7,    0,   14,    7],
       [   0,    8,    0,   19,    9],
       [   0,   12,    0,   26,    8],
       [   0,   15,    0,   31,   10],
       [   0,   17,    0,   38,   15],
       [   0,   16,    0,   47,   17],
       [   0,   21,    0,   53,   23],
       [   0,   29,    0,   63,   28],
       [   0,   36,    0,   91,   39],
       [   0,   51,    0,  102,   47],
       [   0,   57,    0,  127,   58],
       [   0,   64,    0,  149,   74],
       [   0,   83,    0,  189,   92],
       [   0,  106,    0,  245,  111],
       [   0,  128,    0,  319,  129],
       [   0,  151,    0,  381,  158],
       [   0,  189,    0,  462,  190],
       [   0,  236,    0,  555,  217],
       [   0,  283,    0,  697,  276],
       [   0,  364,    0,  893,  351],
       [   0,  436,    0, 1134,  442],
       [   0,  557,    0, 1401,  581],
       [   0,  702,    0, 1831,  740]])

Population regulation

Those populations grew exponentially, if they didn't die out. What if we want a stable population? Let's say that there's "room" for around $K$ individuals, total, and so that each new individual finds "space" to live with probability $1 - N/K$ when there's $N$ adults present. So,

1. Each individual produces a random, Poisson number of offspring, with mean $\lambda$, and then
2. each adult survives with probability $p$ (or else, dies), while
3. each new offspring survives with probability $1-N/K$, if there are $N$ individuals currently.
4. Every individual does this independently of everything else.

This is (a stochastic version of) the logistic model of population growth, with fecundity/recruitment/juvenile mortality regulation.

def logistic_generation(N, lam, p, K):
    """
    N: current number of individuals
    lam: mean number of offspring per capita
    p: probability of survival per capita
    K: carrying capacity (total amount of spaces for individuals)
    Returns: the number of individuals in the next time step
    """
    births = np.random.poisson(lam, N)
    survivors = (np.random.uniform(size=N) < p)
    surviving_births = (np.random.uniform(size=sum(births)) < 1 - N/K)
    return sum(survivors) + sum(surviving_births)

# run five simulations
N = np.full((200, 5), -1, dtype='int')
N[0, :] = 10
for t in range(1, len(N)):
    for k in range(N.shape[1]):
        N[t, k] = logistic_generation(N[t-1, k], lam=2, p=0.5, K=100)


assert(np.min(N) >= 0)

fig, ax = plt.subplots()
ax.plot(N)
ax.set_xlabel("time")
ax.set_ylabel("pop size")
ax.set_title("Logistic growth")

Text(0.5, 1.0, 'Logistic growth')

Exercise: explore parameter space

How do the three parameters, $\lambda$, $p$, and $K$, affect the dynamics? For instance: the stable population size? the rate of growth when small? Try to explain what and why.

Math: it works

Now we'll see how our math-based predictions match the "reality" of our simulations.

Speeding this up with probability

But first, let's speed up the code a bit, so we can simulate larger populations.

Above we did sum(np.random.uniform(0, 1, N) < p) to find the number of survivors: all N individuals flip p-coins to decide if they live or not. Why, that's just the Binomial($N$, $p$) distribution! So, we could have just done np.random.binomial(N, p, 1).

Also, we did sum(np.random.poisson(lam, N)). But, what about additivity of Poissons? If $X_1, \ldots, X_n$ are all Poisson random variables, with means $\lambda_1, \ldots, \lambda_n$ respectively, then $X_1 + \cdots X_n$ is also Poisson, with mean $\lambda_1 + \cdots + \lambda_n$. That means we could have done np.random.poisson(N * lam, 1) instead.

Finally, we will be producing a Poisson number of offspring, then each one flips a coin to decide whether to survive. By the Poisson thinning property, the resulting random number is again Poisson distributed. In other words, if $N$ is Poisson($\lambda$) and $M$ is Binomial($N$, $q$) given $N$, then $M$ is Poisson($\lambda q$).

We'll also vectorize this, to do many simulations at once.

def logistic_gen(N, lam, p, K):
    # N can be a vector of current simulation states
    return (np.random.poisson(lam * N * np.fmax(0, 1 - N / K), len(N)) 
            + np.random.binomial(N, p, len(N)))

logistic_gen(N=np.array([10, 20, 30, 40, 50]), lam=2, p=0.5, K=20)

array([12,  6, 16, 17, 26])

It's equation time

What we're seeing can be thought of as a combination of two things: a deterministic trend, with randomness on top. To describe this, we will work out the mean and variance of the population size at the next time step, as a function of the current population size: $$\begin{aligned} \mathbb{E}[N_{t+1} | N_t = n] &= F(n) \\ \text{Var}[N_{t+1} | N_t = n] &= V(n) . \end{aligned}$$

Suppose that the current population size is $N_t = n$. Given this, what's the mean and variance of $N_{t+1}$? Well, if $S$ is the number of survivors, and $T$ is the number of offspring, then $N_{t+1} = S + T$, and $$\begin{aligned} T &\sim \text{Poisson}(\lambda n (1-n/K)) \\ S &\sim \text{Binomial}(n, p) . \end{aligned}$$ (The distribution for $T$ follows from the additivity and thinning properties of the Poisson distribution, above.) These have mean and variance $$\begin{aligned} \mathbb{E}[T] = \text{Var}[T] = \lambda n (1-n/K) \\ \mathbb{E}[S] = np \qquad \text{Var}[S] = np(1-p) . \end{aligned}$$

Putting this together, we have that $$\begin{aligned} F(n) &= n \left( p + \lambda \left(1 - \frac{n}{K}\right)\right) \\ &= n + n \left( p + \lambda - 1 - \frac{\lambda n}{K}\right) \\ V(n) &= n \left( p (1-p) + \lambda \left(1 - \frac{n}{K}\right) \right) \\ &= F(n) - n p^2 . \end{aligned}$$

The equilibrium

If $F(n_*) = n_*$, then the expected size of a population in the next generation that starts with $n_*$ individuals is unchanged. If not for randomness, a population starting with $n_*$ individuals would stay at that size forever. To find if such a population size exists, and if so, what it is, we solve $$ F(n_*) - n_* = n_* \left( p + \lambda - 1 - \frac{\lambda n_*}{K} \right) = 0 . $$ This has two solutions, zero and $n_* = K (\lambda + p - 1)/\lambda$. If this second solution is positive, then there's a positive equilibrium value.

def logistic_equil(lam, p, K):
    # return predicted equilibrium
    return K * (lam + p - 1) / lam

Time to generalize

And, to save on typing, here's code to "record a bunch of simulations in the columns of an array".

def run_sim(N0, gen_fn, ngens, dtype='int', **kwargs):
    N = np.empty((ngens, len(N0)), dtype=dtype)
    N[0, :] = N0
    for t in range(1, ngens):
        N[t, :] = gen_fn(N[t-1, :], **kwargs)
    return N

Lots of simulations

Ok, finally let's run the simulations, now with our theoretically predicted equilibrium size on top.

N = run_sim(np.repeat(1, 100), gen_fn=logistic_gen, ngens=500, 
            lam=0.5, p=0.81, K=1000)

fig, ax = plt.subplots()
ax.plot(N)
ax.axhline(y=logistic_equil(lam=0.5, p=0.81, K=1000))
ax.set_xlabel("time")
ax.set_ylabel("pop size")
ax.set_title("Logistic growth")

Text(0.5, 1.0, 'Logistic growth')

Noise?

Notice above that the variance of the next time step is proportional to the the number of deaths plus the number of births. This means the standard deviation is proportional to the square root of this. Suppose that in a population of size $N$ there are $R$ births plus deaths. This means that roughly a proportion $R/2N$ of the population is getting replaced each time step. At the next time step, we expect the population to be of size $F(N)$ plus or minus some "noise" of size $\sqrt{R}$. Since $R$ can't be much bigger than $N$, $\sqrt{R}/N$, the size of the noise relative to the population, is of order $1/\sqrt{N}$.

This means that we expect the step-to-step contribution of noise to be smaller in bigger populations.

We can see that in the plot above: the population is roughly of size $N=100$, and the fluctuations about the average are of order 10.

Homework 1

The Ricker model replaces $1 - N/K$ with $\exp(-N/K)$. How similar is it? Simulate from it, solve for the equilibrium value, and check your answer with simulations.

Logistic cobwebs

Now let's look more at the function $n \mapsto F(n)$. We can do this the "math-free" way by just plotting the population size after one time step of the simulation starting from a grid of possible initial sizes. The diagonal red line is where $F(n) = n$, i.e.: equilibrium!

X = np.arange(250)
Y = logistic_gen(X, lam=0.5, p=0.95, K=200)

plt.plot(X, Y)
plt.plot([0.0, 250], [0.0, 250])
plt.xlabel("N(t)")
plt.ylabel("N(t+1)")
plt.show()

That's hard to see! Let's look at $n \mapsto (F(n) - n)$ instead. The horizontal line is where $F(n) - n = 0$, again: equilibrium!

plt.plot(X, Y-X)
plt.axhline(0.0)
plt.axvline(x=logistic_equil(lam=0.5, p=0.95, K=200), color='r')
plt.xlabel("N(t)")
plt.ylabel("N(t+1) - N(t)")
plt.show()

New parameters

Now let's spice things up by changing the parameters to $$\begin{aligned} \lambda &= 5 \\ p &= 0.25 \\ K &= 200 \end{aligned}$$

N = run_sim(np.repeat(170, 20), gen_fn=logistic_gen, ngens=10,
            lam=5.0, p=0.25, K=200)
equil = logistic_equil(lam=5.0, p=0.25, K=200)
    
plt.plot(N, label='logistic growth')
plt.axhline(y=equil)
plt.xlabel("time")
plt.ylabel("pop size")
plt.show()

Wow, what happened? The population is oscillating. Let's look at $F(n)$ to understand.

X = np.fromiter(range(300), dtype='int')
Y = logistic_gen(X, lam=5., p=0.25, K=200)

fig, ax = plt.subplots()
ax.plot(X, Y)
ax.plot([0.0, 300], [0.0, 300])
ax.axhline(y=equil, color='m', linestyle="dotted")
ax.axvline(x=equil, color='m', linestyle="dotted")
ax.set_xlabel("N(t)")
ax.set_ylabel("N(t+1)")
plt.show()

Question: What's that wierd tail?

The cobweb plot

Using a plot of the expected value curve $f(n) = \mathbb{E}[N_{t+1} | N_t = n]$, you can trace the expected dynamics as $$ n \mapsto f(n) \mapsto f(f(n)) \mapsto f(f(f(n))) \mapsto \cdots $$ by "bouncing off the diagonal".

It's easier to explain with a picture.

def cobweb_plot(N, ax):
    for k in range(N.shape[1]):
        X = np.repeat(N[:, k], 2)[:-1]
        Y = np.repeat(N[:, k], 2)[1:]
        ax.plot(X, Y)

# NOTE: run this again to get more cobwebs!
N = run_sim(np.array([equil]), gen_fn=logistic_gen, ngens=20,
            lam=5.0, p=0.25, K=200)     
print("Population size:", N.T)

cobweb_plot(N[:, :1], ax)
fig

Population size: [[170 181 149 220  50 209  48 196  66 219  62 219  60 225  56 235  67 258
   61 229]]

Exercise: dizzy spiders

Fact: If the slope of $F()$ (the blue curve above) at the equilibrium (the magenta dotted lines) is steeper than $-1$, then the equilibrium is unstable (and so the population will tend to oscillate around it).

Experiment to convince yourself of this, draw pictures on the whiteboard with your neighbor, then explain why this happens.

What about noise?

Above we saw that bigger populations had less noise (in this model anyhow). This says that bigger populations should behave more like solutions to the differential equation. Let's see how that holds up in our simulations.

for K in [200, 2000, 20000]:
    print(K)
    N = run_sim(np.array([0.6, 0.8, 1.0, 1.2, 1.4]) * logistic_equil(lam=0.1, p=0.95, K=K), 
                gen_fn=logistic_gen, ngens=120, 
                lam=0.1, p=0.95, K=K)
    plt.plot(N)
    plt.axhline(y=logistic_equil(lam=0.1, p=0.95, K=K))
    plt.xlabel("time")
    plt.ylabel("pop size")
    plt.title("K={}".format(K))
    plt.show()

200

2000

20000

Oscillations, and Chaos

It would be a crime to talk about the logistic equation and not see period doubling, and chaos. A nice review of complicated behavior in simple models is May 1976. In this section we'll just work with the determinstic model (remember that $r = \lambda + p -1$ and we'll switch out $K$ for $n_* = K r / \lambda$, to be consistent with the literature).

def logistic_step(N, r, C):
    return N + r * N * (1 - N / C)

N = run_sim([80], gen_fn=logistic_step, dtype='float',
            ngens=50, r=1.1, C=100)

plt.plot(range(N.shape[0]), N, linestyle="dotted")
plt.scatter(range(N.shape[0]), N)
plt.xlabel("time")
plt.ylabel("pop size")
plt.show()

To see if the model approaches a stable, oscillating equilibrium, we can (a) iterate the model for a while, and then see if the values repeat. For instance, at $r=2.1$ there's a period-2 loop:

N = run_sim([60], gen_fn=logistic_step, dtype='float',
            ngens=50, r=2.1, C=100)

plt.plot(range(N.shape[0]), N, linestyle="dotted")
plt.scatter(range(N.shape[0]), N)
plt.xlabel("time")
plt.ylabel("pop size")
plt.show()

Here's another example (at $r=2.5$, with period 4):

N = run_sim([10], gen_fn=logistic_step, dtype='float',
            ngens=100, r=2.5, C=100)

plt.plot(range(N.shape[0]), N, linestyle="dotted")
plt.scatter(range(N.shape[0]), N)
plt.xlabel("time")
plt.ylabel("pop size")
plt.show()

... and, what's going on here???

N = run_sim([10], gen_fn=logistic_step, dtype='float',
            ngens=100, r=2.9, C=100)

plt.plot(range(N.shape[0]), N, linestyle="dotted")
plt.scatter(range(N.shape[0]), N)
plt.xlabel("time")
plt.ylabel("pop size")
plt.show()

Ok, now let's see how that behaves as we change $r$. To save space, we'll (a) cut out the first bunch of steps of each simulation, and (b) plot the values of $N_t$ against their respective values of $r$ (so the first plot above would be compressed to just two points, at about 82 and 111, above $r=2.1$.

rvals = np.linspace(1.5, 3.0, num=400)
N = run_sim(np.repeat(10.0, len(rvals)), gen_fn=logistic_step,
            dtype='float', ngens=200, r=rvals, C=100)

plot_gens = 20
plt.scatter(np.row_stack([rvals]*plot_gens), N[-plot_gens:,:],
            marker='.')
plt.show()

The first "period doubling", where the single stable cycle becomes unstable, occurs at $r=\sqrt{6}$. The period doubles again, followed soon by the onset of chaos.

Homework 2

Here is another "logistic" model. Suppose that everyone dies every time step, and that large populations deplete the amount of resources left for the next generation. The net effect is that the mean fecundity of an individual in a population of size $N$ is $1 + r - N/\beta$, so that $1+r$ is the mean fecundity in the presence of lots of resources, but that mean fecundity decreases by $1/\beta$ for each additional neighbor. Implement a simulation of this process (using Poisson offspring and Poisson additivity), and plot how it evolves. Write down and solve the equations to find its equilibrium in terms of $r$ and $\beta$.

Bonus: We can't see chaos in our simulations above because of the max(0.0, 1-N/K) term. We can see it in this model. Find parameter values with (a) a bistable cycle and (b) chaos.